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In the stable marriage problem A*' men and A'^ women have 
to be matched by pairs under the constraint that the result- 
ing matching is stable. We study the statistical properties of 
stable matchings in the large A'^ limit using both numerical 
and analytical methods. Generalizations of the model includ- 
ing singles and unequal numbers of men and women are also 
investigated. 
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I. INTRODUCTION 

The study of complex systems that are built up of 
many individuals interacting according to relatively sim- 
ple laws has recently attracted a lot of interest among 
physicists. Concepts and methods developed in statisti- 
cal physics have successfully been applied to describe col- 
lective and nonlinear phenomena arising in such systems. 
Studies on traffic flow, stock markets, voting, flocking 
birds, evolution etc. are only a few examples to mention 
here. 

In the stable marriage problem the elements of two sets 
have to be matched by pairs. As possible applications one 
could think of job seekers and employers, lodgers and 
landlords, or simply men and women who want to get 
married. For the sake of clearness and simplicity we will 
exclusively refer to the paradigm of marriage in this pa- 
per. The purpose of most studies of the stable marriage 
problem has so far been to unravel the underlying math- 
ematical structure and to develop new or improve exist- 
ing algorithms for finding stable marriages Only 
recently it has been recognized that the stable marriage 
problem has many features in common with classical dis- 
ordered systems in statistical physics like spin glasses, in 
particular the element of frustration and the existence 
of many competing states [||. Applying methods devel- 
oped in statistical mechanics, Nieuwenhuizen [Q solved 
several variants of the stable marriage problem, including 
bachelors and polygamy. In his approach, Nieuwenhuizen 
focussed on the properties of globally optimal matchings 
which are advantageous for the society as a whole, but 
not necessarily for all individuals. Here we consider the 
opposite limit where all individuals pursue only their own 
egoistic objectives regardless of the consequences for oth- 
ers. Under these conditions, only matchings that are 



stable with respect to the actions of all individuals are 
relevant. In the framework of game theory this kind of 
stability is termed Nash equilibrium. 

Our paper is organized as follows. In Section 2 we 
introduce the stable marriage model and set up the no- 
tation. The Gale-Shapley algorithm and the properties 
of the stable matchings obtained with it are presented in 
Section 3. In Section 4 we analyze the set of all stable 
matchings. Generalizations of the model, introducing an 
acceptance threshold and unequal number of men and 
women are discussed in Sections 5 and 6, respectively. 
Finally, we summarize our results in Section 7. 



II. THE STABLE MARRIAGE MODEL 

We consider a society consisting of iV men and N 
women who have to be matched by pairs. The cost for 
man i to be married with woman j is x{i,j) and the 
cost for woman j to be married with man i is y{j,i)- 
In order to investigate statistical properties of the model 
we assume that all cost functions are independent ran- 
dom variables uniformly distributed between and 1 0. 
Since all individuals simultaneously attempt to optimize 
their own benefit by being married with a partner of low- 
est possible cost, conflicting wishes arise, and in general 
it is impossible to find a matching where everybody is 
perfectly satisfied. Let us assume that a matching is 
formed where man i is married with woman j. We define 
Xi = x{i,j) and ijj = y{j,i), i.e. Xi is the cost for man i 
in this particular matching and yj the corresponding cost 
for woman j. Instead of cost we will frequently use the 
expression energy in order to emphasize the analogy with 
physical systems where the energy is the quantity to be 
minimized in the groundstate. For later convenience, we 
introduce the total energies of men and women by 



X 



N 



Y 



N 



(1) 



In search of the rules according to which marriages can 
be formed or broken, essentially two points of view have 



^In Ref. Q a different probability distribution has been con- 
sidered, namely p{x) oc x"^ exp(— a;). The case r = corre- 
sponds essentially to our choice. 
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emerged: In the first one, the goal is to find the glob- 
ally best matching, i.e. the one where the total energy 
U ^ X + Y is minimal. This problem has recently been 
studied by Nieuwenhuizcn |^ using techniques borrowed 
from statistical physics of disordered systems, such as 
the replica trick. The second point of view, that we will 
adopt, emphasizes the role of individuals as decision- 
makers. As long as there is no supervising body like 
church or government that dictates who has to be mar- 
ried with whom no individual is obliged to stay in an 
unsatisfactory relationship if he or she can find a better 
one. This leads us to the notion of stability. A matching 
is called stable if there is no man and no woman who pre- 
fer each other to their actual partners. In other words, a 
matching is unstable, if there exists a man i and a woman 
J, who are not married with each other, such that 



> x{ij) and > y{j,i) 



(2) 



The particular significance of stable matchings lies in the 
fact that they cannot be broken, neither by the action of 
a single pair nor by any coalition of men and women. In 
contrast to the globally best matching which can only be 
found with the help of computers for very small numbers 
of men and women, there exist powerful algorithms to 
search for stable matchings, even for systems of several 
thousands of individuals. The simplest of these algo- 
rithms, introduced by Gale and Shapley some time ago 
M , will be presented in the next section. 



III. THE GALE-SHAPLEY ALGORITHM 

In the man-oriented Gale-Shapley (GS) algorithm only 
men make proposals while women accept or reject. In the 
beginning every man makes up a list where he places the 
women according to his preferences. The smaller the cost 
x{i,j) the higher the position of woman j on the list of 
man i. Then a series of proposals starts. Whenever a 
man is not married he makes a proposal to the woman 
with the highest rank among those to whom he has not 
yet proposed. If the woman is not engaged she accepts 
the proposal and they get married. If the woman is al- 
ready married, she accepts if she prefers the suitor to her 
husband. In this case she divorces and gets married with 
the proposer. The algorithm stops when the last woman 
has received an offer and gets married. The number of 
proposals r that a man has to make on the average during 
the execution of the GS algorithm has been calculated in 
§1. The result is 



r~\ogN + G 



(3) 



where C = 0.5772... is Euler's constant and corrections 
are negative and of order (logA^)^/A^. Amazingly, the 
knowledge of r is already sufficient to calculate the total 
energy of men (A) and women (Y) in the GS matching. 
With every proposal the energy of one man increases by 



~ 1/N since he has to move forward one step on his 
preference list. Consequently, the total energy of men is 
on the average the same as the number of proposals that 
every man makes, i.e. X = r. 

Each woman also receives on the average r propos- 
als. Keeping only the best offer i.e. the smallest out of 
r random numbers between and 1 yields a best value 
of y = l/(r -|- 1). Taking into account that the num- 
ber of proposals that a woman receives is not fixed, but 
distributed according to a binomial distribution yields 
y ~ {1 — exp(— r))/r w l/r (see Appendix B). Using Eq. 

we obtain the following estimates for the energies of 
men and women in the man-oriented GS matching 



X~loeN + C Y 



N 



logA^-l-C 
which implies the relation XY^N. 



(4) 
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FIG. 1. Average energy of men and women in the 
man-oriented Gale-Shapley algorithm as a function of N. 

Fig. ^ shows the average energies X and Y of men and 
women obtained with the man-oriented GS algorithm as 
function of N. Each data point represents an average 
over several thousand realizations. A logarithmic scale is 
used for the energy axis in order to highlight the good 
agreement with the analytical result of Eq.(^ even for 
small values of N. 

Of course, everything that has been said about the 
man-oriented GS algorithm is also valid in the woman- 
oriented version where the roles of men and women are 
interchanged, i.e. women propose and men judge. It has 
been shown |^ that the stable matching obtained from 
the man-oriented execution of the GS algorithm is man- 
optimal in the sense that no man can have lower energy 
in any other stable matching. Similarly, it is the worst 
possible stable matching for women. Therefore the ener- 
gies obtained by performing the GS algorithm provide an 
upper and a lower bound for the possible energy range 
of all stable matchings. Furthermore, if the man and 
the woman-optimal matching happen to be identical we 
know for sure that in this case there exists only one single 
stable matching. 
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IV. THE SET OF ALL STABLE MATCHINGS 

If N is large, there are generally many other sta- 
ble matchings besides the man and the woman-optimal 
matchings that can be obtained with the GS algorithm. 
Although the maximum number of stable matchings in a 
system of size N is not known, there exist lower bounds 
1^; e.g. for = 32 the maximum number of stable 
matchings is larger than 10^^ (!). 

In order to calculate the average number of stable 
matchings we consider an arbitrary matching of N men 
and N women. This matching is unstable, if there exists 
a man i and a woman j who are not married with each 
other, and whose mutual costs x{i, j) and y{j, i) are both 
lower than the ones in the existing matching, Xi and yj, 
respectively. The probability that man i and woman j 
prefer to stay with their respective partners is therefore 
Pij = 1 — Xi yj and the whole matching is stable with 
probability 



P 



d^x 



11(1 



Xi yj) 



(5) 



Although this 2A^-fold integral, originally derived by 
Knuth M, looks quite simple, it cannot be done exactly 
unless N is very small. Using rather involved proba- 
bilistic arguments, Pittel has derived an asymptotic 
formula, 



logA^ 
eT{N) 



(6) 



valid in the limit N oo. However, comparing Eq. 

with numerical simulations we found rather serious 
discrepancies (see Fig. 2). In order to remove these dis- 
crepancies and to provide a more transparent approach 
we replace the factors 1 — Xi yj in Eq. (||) by exp(— a;^ y^), 
which is justified since the dominant contribution to the 
integral comes from the region where the products Xiyj 
are small. Dropping the constraint i 7^ j we arrive at the 
approximation 



P = 



d^y exp(-Xy) 



(7) 



where the variables X = X^ili ^^'^ ^ — '^f=i Vj 
the energies of men and women, respectively. Eq. (|^) 
can also be viewed as the partition sum of two species of 
particles confined to the interval [0, 1] whose interaction 
is the product of their center of mass coordinates. It is 
convenient to replace the integral over the variables Xi 
and yj by an integral over X and Y 



P = 



dX dY p{X)p{Y) exp{-XY) (8) 



'0 Jo 

where the probability distribution of X is asymptotic to 



X 



N-1 



r{N) 



l-exp(- — ) 



N 



(9) 



for small X (see Appendix A). Inserting Eq. (^ in Eq. 
(^) and introducing t = XY as new integration variable 
we obtain 



p = r{Ny 



dxx- 



N X 
log N 



dt t^-^ exp(-i) 



T{N)-^ i dXX-^ 

JlogN 

r{N)-^ (logA^- 21og(logA^)) 



(10) 



N 



where we have assumed that the factor (l — exp(— ^)) 
imposes a cut-off ~ N/ log N on the X-integration and 
the same holds for Y. The integral over t yields r(AJ, 
provided that X > \ogN. Comparison with Eq. (n) 
shows that the ratio P/P approaches e = exp(l) in the 
limit A^ — > 00 which we also checked numerically. Since 
there are in total A^! = r(A^ -|- 1) matchings the average 
number of stable matchings is given by 



N 

5 ~ — (log A - 2 log(log A^)) 
e 



(11) 




FIG. 2. Average number of stable matchings S as function 
of TV. Numerical simulations (circles) in comparison with the 
analytic expression of ref. [6] (solid curve) and our result of 
Eq.(|ll|) (dashed curve). 

S is plotted as function of A in Fig. |^. The circles 
are data obtained from numerical simulations using the 
algorithm of Gusfield and Irving |^ that allows the deter- 
mination of all stable matchings. Each data point repre- 
sents an average over 200 random realizations of the pref- 
erence lists. The dashed curve is the result of Eq. ( |Tl| ) 
and the solid curve is the asymptotic formula derived by 
Pittel H]. Note the large error made by neglecting the 
log(logA^) term. 

In the following we investigate how the energies of men 
and women are correlated in stable matchings. From Eq. 
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(|) we see that p{X, Y) = p{X)p{Y) exp(-Xr) can be 
considered as the probabiHty density of finding a sta- 
ble matching with energies X and Y. Due to the factor 
(l — exp(— ^)) that appears in p{X) the occurrence 
of stable matchings with X > N/ log N is very unlikely, 
and by symmetry, the same cutoff exists for Y. This ob- 
servation is in accordance with our analysis of the GS 
algorithm, where we found that the energy of women in 
the man-optimal stable matching (i.e. the worst possi- 
ble for women) is ~ A^/ log TV, and vice versa. In the 
region log N < X,Y < N / log N the probability distri- 
bution simplifies to p{X,Y) cx (Xy)^^i exp(-Xy). It 
depends only on the product of men's and women's ener- 
gies and is sharply peaked around the curve XY = N for 
large values of N . Mean- value and standard deviation of 
the product XY are given by 



{XY) = TV, (t{XY) 



N 



(12) 



where (. . .) denotes the average calculated with p{X, Y). 
At the symmetric point X = Y = y/N the width of the 
distribution is 0(1) which means that relative fluctua- 
tions go to zero for — > oo. 



^ 0.4 
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FIG. 3. Probability distribution of log 5* for = 50. Nu- 
merical simulations (circles) fitted by a Gaussian. 




FIG. 4. Standard deviation a of log S as function of A'^. 
The numerical results (circles) seem to converge to the solid 
line (Too « 0.48 for A' ^ oo. 

Finally, we investigate the probability distribution of 
the number of stable matchings using numerical simula- 
tions. While S has a very asymmetric distribution which 
decays quite slowly for large values of S, the distribution 
of log S can be nicely fitted with a Gaussian, as shown in 
Fig. ^ for A^ = 50. Remarkably, the width of the distri- 
bution, i.e. the standard deviation a{\og S) approaches a 
constant value aoo ~ 0.48 for large A^ as shown in Fig. |^. 
We have no argument to explain this universal behavior. 



V. ACCEPTANCE THRESHOLD 

Up to now we have assumed that every person has to 
get married. In the following we soften this condition by 
introducing an acceptance threshold A < 1 which is the 
energy of individuals who remain single. In principle, 
A could be chosen at random for all individuals, or at 
least be different for men and women, but for the sake of 
simplicity we consider only the case where A is the same 
for all men and all women. 

Let us consider a situation where A is sufficiently large 
that singles do not occur. The GS algorithm is readily 
generalized to the case with acceptance threshold. In 
the man-oriented version the only difference is that now 
women - even when unmarried - refuse proposals from 
men whose energy is higher than A. This increases the 
average number of proposals that are necessary to arrive 
at a stable matching by a factor of 1/A. The generaliza- 
tion of Eq. (0) to the case A 7^ 1 reads 



X{A) 



\ogN + C 
A 



y(A) ~ A 



N 



logA^-t-C 



Both energies coincide for 

logN + C 



N 



(13) 



(14) 



and we conclude that for A < Ac there exist only few 
or one single stable matching. This fact can be used to 
generate approximately sex-fair stable matchings - which 
is difficult to achieve by other means - using the GS 
algorithm and setting A = Ac. 

In the limit when A is sufficiently small it becomes 
favorable for a certain fraction of men and women to re- 
main single. In this case it is hard to calculate the average 
number of proposals r explicitly since the series of pro- 
posals is not only terminated when all men are married 
but also if no acceptable women are available. On the 
other hand the energies of men and women coincide in 
this regime since there exists only one stable matching. 
Let us assume that a woman has received r proposals 
out of which she keeps only the best one. If none of the 
proposals are better than the threshold A she remains 
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single at the cost of A. Analyzing the probability dis- 
tribution of the best offer (see Appendix B) yields the 
average energy of women 



y = I (1 - exp(-XA)) 

From the condition X = Y we obtain 
^N{l-exp{-XA)) 



(15) 



(16) 



In general this equation can only be solved numerically. 
The limiting cases are 



X(A) 




(17) 



which means that for sufficiently small A it is advan- 
tageous for almost everybody to remain single (at the 
energy A) while for A > 1/ \/N the energies in sex- fair 
matchings are X = Y = \/N, in agreement with our 
previous considerations. 




A 

FIG. 5. Energy of men as function of A. Open circles: 
solution of lowest energy (man-optimal GS), filled circles: so- 
lution of highest energy (woman-optimal GS). The solid and 
the dashed curves are the analytical results of Eq. ( |l^ ) and 
Eq. (|l^), valid for A > Ac and A < Ac, respectively. 

In Fig. ^ we show the energy of men and women in 
the man-optimal stable matching as function of A for 
N = 200. Each data point represents the average over 
1000 executions of the GS algorithm. The agreement be- 
tween numerical and analytical results is very good with 
the exception of the cross-over region around Ac where 
neither of the conditions that were made to treat the 
problem (absence of singles in the case A > Ac, exis- 
tence of only one stable matching in the case A < Ac) is 
strictly fulfilled. 

In the following we investigate how the number of sin- 
gles depends on the acceptance threshold A. The prob- 
ability that a woman who receives r proposals rejects all 
of them and remains single is ps = (1 — A)''. Taking 
into account that the number of proposals addressed to 



a particular woman is distributed according to a binomial 
distribution the average number of singles is related 
to the energy X via 



Ns{A) = N exp(-A:A) 



(18) 



(see Appendix C) and can be calculated explicitly from 
the solution of Eq. (|l6|) . The limiting cases are 



J N exp(-A2iV) 
^«(^) - { N exp(-AViV) 



for 
for 



(19) 



When A is decreased, the first single appears when 
Ns{A) = 1. This is the case for A = log N/VN ~ Ac in 
accordance with our assumption that there are no singles 
for A > Ac. 

Finally, we generalize our results concerning the aver- 
age number of stable matchings to the case with thresh- 
old A < 1. The only difference compared to the calcula- 
tion for A = 1 is that the probability distribution of X 
(see Appendix A) has to be replaced by 



piX) 



X 



N-1 



TiN) 



1 — cxp(- 



iVA, 



N 



(20) 



which takes into account that all energies Xi have to be 
smaller than A. Making the same approximations as be- 
fore we obtain for the average number of stable matchings 



^(A) 



N 



(logiV-21og(logiV) + 21ogA) (21) 



S'(A) goes to zero at A = log N/y/N ~ Ac in agree- 
ment with the fact that man and woman-optimal stable 
matchings coincide for A < Ac. 



VI. UNEQUAL NUMBER OF MEN AND WOMEN 

Let us now consider the more general case where the 
number of men and women is not the same. To be spe- 
cific, we assume that there are A^-l- 1 men and A'' women. 
It is obvious that in this case at least one man has to 
remain single. As before, the GS algorithm can be used 
to find the man-optimal and the woman-optimal stable 
matching, respectively, depending on who proposes and 
who judges. There is now a fundamental asymmetry be- 
tween the man and the woman-oriented version: when 
women propose the algorithm stops as usual when the 
last woman has found a husband; in the man-oriented 
version, however, the series of proposals stops only when 
the first man has reached the bottom of his list, i.e. when 
he has been rejected by all women and must stay sin- 
gle. Therefore, the members of the minority group (the 
women in our case) are substantially favored. 
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FIG. 6. Energy of men (circles) and women (triangles) as 
function of A for 200 women and 201 men. Full and open 
symbols denote the best and the worst possible solution, re- 
spectively. The solid and the dashed curves are the analytical 
results discussed in the text. 

This is illustrated in Fig. ^ where the energy of men 
(circles) and women (triangles) is plotted versus the ac- 
ceptance threshold A. The data points represent aver- 
ages of 1000 samples using the GS algorithm for 201 men 
and 200 women. Best (open symbols) and worst solu- 
tions (filled symbols) are much closer in energy now than 
in the symmetric case with equal number of men and 
women (see Fig. ^ for comparison) indicating that the to- 
tal number of stable matchings is reduced. The energies 
in all other stable matchings lie in the two narrow corri- 
dors between the best and the worst solution. The dashed 
curves are X(A) = X(l)Aandy(A) = r(l)/A obtained 
from the generalization of Eq. (|l^) to the present case 
while the solid curve is the solution of Eq. ( p!6| ) valid for 
A < Ac. 



VII. CONCLUSIONS 

We have investigated the statistical properties of sta- 
ble marriages using both numerical and analytical tools. 
A detailed analysis of the GS algorithm is provided. This 
algorithm yields two extreme stable matchings and is 
thus particularly useful to characterize the whole set of 
stable matchings. The average number of stable match- 
ings is calculated from an integral representation that can 
be viewed as the partition sum of two species of particles 
living on the unit interval, with multiplicative interaction 
of their center of mass coordinates. Our analytical result 
agrees very well with numerical simulations, in contrast 
to a previously obtained, less accurate formula. In addi- 
tion, the fundamental relation XY = N between the en- 
ergy of men and women in stable matchings follows quite 
naturally from our approach. Numerical simulations in- 
dicate that the probability distribution of the logarithm 
of the number of stable matchings is Gaussian with a 
universal width CToo ~ 0.48 for N — > oo. 



When an acceptance threshold A is introduced, the 
number of stable matchings shrinks and singles occur. 
Three different regimes can be identified: For A > Ac ~ 
log N/^/N there are many stable matchings whose en- 
ergies extend over a wide range and no singles occur. 
For 1 / x/iV < A < Ac there are very few or only one 
stable matching. In this regime the energy of men and 
women is nearly independent of A although the num- 
ber of singles increases exponentially with decreasing A. 
For A < there are many singles and only a small 

fraction of married individuals. When the number of 
men and women is not the same, the number of stable 
matchings is also substantially reduced. Furthermore, 
the members of the minority group are strongly favored 
compared to the majority group, even if there is only one 
additional member. 

The general ideas discussed in this paper could become 
important for investigations on the labour market since 
nowadays information about forthcoming job opportuni- 
ties, on the one hand, and qualifications of applicants, 
on the other hand, can be made available to everybody 
via the internet. Preliminary results on applications of 
the stable marriage problem in the context of the labour 
market can be found in |^] . Further work on this subject 
is in progress. 
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APPENDIX A: CALCULATION OF p{X) 

In this appendix wc calculate the distribution function 
p{X) of the variable X — Xi which is formally given 

by 



p{X)= ['d'^xSiX-Y, 
Jo 



Xi) 



(Al) 



According to the central limit theorem p{X) converges 
to a normal distribution 



piX) = 



1 



V27: 



exp(~^^f#^) (A2) 



with X = N/2 and cr^ = N/12 for large N. This is 
however only true for the central part of the distribution 
whereas the tails of p{X) look very different. Let us 
introduce new variables qi such that Xi = qt — qi-i for i = 
1,N and qo = 0,qN = X . Disregarding the constraint 
Xi < 1 for the moment we obtain 



piX) 



1 



X 



dqi 



X 



dgjv-i 



X 



T{N) 



(A3) 
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where the factor r(A'') occurs due to the fact that the 
qi have to be arranged in ascending order to satisfy 
the condition Xi > 0. The condition Xi < 1 can 
now be to implemented as foUows: Considering the qi 
as independent random variables uniformly distributed 
in the interval [0, X] we obtain a Poisson distribution 
p{x) — Aexp (—Ax) for the distance x between two con- 
secutive numbers, where A = N/X. The probability that 
all < 1 is therefore given by 



P(aU Xi < 1) 



1 N ^ 

dx p{x) 







1 — exp (- 



iV, 
X 



N 



(A4) 



neglecting the weak correlations among the Xi d ue t o the 
fact that their sum is fixed. Combining Eqs. (A3) and 



(A4) we obtain the final result 



p{X) 



X 



N-1 



valid for X < N/\ogN. 



exp(-^. 



N 



(A5) 



APPENDIX B: AVERAGE ENERGY OF WOMEN 

Let us assume that men propose in the GS algorithm 
and that a woman receives r proposals out of which she 
keeps only the best one. If none of the offers is better 
than A she remains single at the energy A. The average 
energy thus obtained is 



1 



V 



1 



(l-(l-A) 



r+l\ 



(Bl) 



We have to take into account that the number r of pro- 
posals that a woman receives is not fixed but distributed 
according to the binomial distribution 



h{T) 



R 



R-r 



(B2) 



where R is the total number of proposals that men make 
during the execution of the GS algorithm and p ~ is 
the probability that such an offer is addressed to a given 
woman. Averaging Eq. (Bl) with respect to h{r) yields 



1 



R 

E 

N 



R 



P^{l-p) 



R-r 



1 



and similarly 

(1- A)'-+i 
r + 1 



R+1 



N 
R+1 



1-(1- 



\R+1 



N' 



(B3) 



A, 



\ at' ^ AT' 



R+1 



N 



(B4) 



Combining Eqs. (Bl), (B3) and (B4) we obtain 



(y) = ^ f 1 - (1 - ^r^' 



N f. / -Ran 

^(l-exp(--A) 
(B5) 



The average energy of women Y — N{y) is therefore 
related to the energy of men X = R/N via 

N 

r = - (1 - exp i-XA)) (B6) 



APPENDIX C: NUMBER OF SINGLES 

The probability for a woman to remain single after 
receiving r proposals is ps = (1 — A)'. Averaging with 
respect to the binomial distribution (B2) yields 



r=0 ^ ^ 



P) 



R-ri 



(CI) 



where p = 1/N and R is the total number of proposals 
made. Using X = R/N we obtain 



{ps) = exp(-XA) 



(C2) 
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